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We experiment with modifications of the BSSN form of the Einstein field equations (a refor- 
mulation of the ADM equations) and demonstrate how these modifications affect the stability of 
numerical black hole evolution calculations. We use excision to evolve both non-rotating and rotating 
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and obtain accurate and stable simulations for specific angular momenta J/M of up to about 0.9A/. 
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I. INTRODUCTION 

Binary black holes are among the most promising 
sources for the gravitational wave laser interferometers 
currently under development, including LIGO, VIRGO, 
GEO, TAMA and LISA. The identification and interpre- 
tation of possible signals requires theoretically predicted 
gravitational wave templates. For the late epoch of the 
binary inspiral, numerical relativity is the most promis- 
ing tool for the computation of such templates. 

The numerical simulation even of single black holes 
has encountered numerous difficulties, which presumably 
arise from the complexity of Einstein's equations, the ex- 
istence of a singularity inside the black hole, and the 
gauge (or coordinate) freedom inherent in general rela- 
tivity. Some recent developments, however, have led to 
significant and very promising advances. 

Traditionally, most numerical relativity simulations 
were based on the 3 + 1 decomposition of Arnowitt, Deser 
and Misner (ADM, Q), which has been shown to de- 
velop instabilities often H Q- To avoid these instabil- 
ities, a number of hyperbolic formulations of Einstein's 
equations have been developed (see, e.g., 0, 3 as well 
as [j| and references therein) . Alternatively, Shibata and 
Nakamura and Baumgarte and Shapiro Q introduced 
a modification of the original ADM equations that in- 
volves a conformal-traceless decomposition and the in- 
troduction of a new auxiliary variable (see Section ITTlbe- 
low) . This formulation, now commonly referred to as the 
BSSN formulation, has led to significant improvements 
over the original ADM equations, and has been widely 
adopted (e.g. 0, H 13, E3, 111 El ) 

Singularities inside black holes were traditionally 
avoided by using "singularity avoiding" coordinate condi- 
tions, including maximal or polar slicing [l3l liH ll5l lTr ^ | . 
Typically, these conditions lead to grid pathologies that 
cause codes to crash after relatively short times. An al- 
ternative strategy takes advantage of the fact that the 
black hole exterior is causally disconnected from the in- 
terior, so that a region of the black hole interior in- 
cluding the singularity can be excised from the com- 
putational grid. These "singularity excision" techniques 



E1II1 EM [Hill have led to large improvements in 
the simulation of single black holes, and are a promising 
tool for binary black hole evolutions (see for prelim- 
inary results). 

The application of singularity excision requires a co- 
ordinate system that is regular across black hole hori- 
zons, allowing smooth horizon penetration. For single 
black holes, one such coordinate system is the Kerr-Schild 
form, which can be used to represent both Schwarzschild 
and Kerr black holes (and which is also invariant under 
boosts). Binary black hole initial data based on Kerr- 
Schild coordinates can be constructed by solving the con- 
straint equations of general relativity for "corrections" 
arising from su perpos itions of two boosted Kerr-Schild 
black holes |H 111 |H |28t 

Alcubierre and Briigmann |29| recently combined the 
BSSN formalism with a particularly simple singularity 
excision method to evolve single black holes in Kerr- 
Schild coordinates. Restricting the evolution to octant 
symmetry, all fields settle down to equilibrium (of the 
finite-difference equations), and no instabilities are en- 
countered. If, however, the symmetry assumption is re- 
laxed, instabilities develop and the code crashes after a 
few hundred M. Similar findings were reported in |30| . 
where a completely independent formulation and imple- 
mentation was adopted. Improvements over these results 
were discussed in [31], |32|, |33| , but even these do not com- 
pletely eliminate the instabilities for evolutions without 
symmetry assumptions. 

Recent results suggest that adding constraints to the 
evolution equations affects the numerical stability of the 
system (e.g. |H IM IM El E3; see jH for an illus- 
tration in electrodynamics). In this paper we experi- 
ment with adding the new constraints that appear in the 
BSSN formulation to the evolution equation of the new 
auxiliary functions, and, following [33, the Hamiltonian 
constraint to the evolution equation for the spatial met- 
ric. We also experiment with schemes for imposing alge- 
braic constraints on the conformally related metric and 
extrinsic curvature, as well as with different shapes for 
the excised region inside the black hole. 

With these modifications we obtain evolutions of single 
black holes that last over several thousand M, indepen- 
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dent of any symmetry assumptions, without encounter- 
ing any evidence of a growing instability. These findings 
hold both for static black holes and for rotating black 
holes with specific angular momentum of up to 0.9M. 

The paper is organized as follows: We summarize the 
BSSN formulation in Sec. [H] and black hole spacetimes 
in Kerr-Schild coordinates in Section ITTT1 Our modifica- 
tions of the BSSN scheme are described in Sec. IIVI In 
Section IVT1 we present results of our simulations for both 
static and rotating BHs. We summarize and discuss the 
implications of our findings in Sec. IVIII We also include 
an Appendix that explains our evaluation of the ADM 
mass and angular momentum. Throughout the paper we 
adopt geometrized units with G = c= 1 . 



with the conformal factor chosen so that the determinant 
7 of jij is unity 



(9) 



The traceless part of the extrinsic curvature K~ij , defined 
by 



(10) 



where K = r y v is the trace of the extrinsic curvature, 
is conformally decomposed according to 



(11) 



II. THE BSSN FORMULATION 



The conformal connection functions T l , initially defined 
as 



We write the metric in the ADM form 

ds 2 = -a 2 dt 2 + jijidx* + (3 l dt){dx j + p j dt), (1) 

where a is the lapse function, (3 l is the shift vector, and 
7y is the spatial metric. Throughout this paper, Latin 
indices are spatial indices and run from 1 to 3, whereas 
Greek indices are space-time indices and run from to 3. 

The Einstein equations can then be decomposed into 
the Hamiltonian constraint Ti. and the momentum con- 
straints Mi 



H = R— + K 2 = 0, 

Mi = DjK\ — DiK = 0, 

and the evolution equations 



dt 13 31 



(2) 
(3) 



(4) 



-K^ = -DiDja + a(Rij - 2K lt K l 3 + KK t] ). (5) 



Here we have assumed vacuum T a p = and have used 
d d 



dt dt 



(G) 



where Cp is the Lie derivative with respect to [3 l . Di is 
the covariant derivative associated with 7^ , Rij is the 
three-dimensional Ricci tensor 



r> kt i \ 

+7 (T m ieT m kj — T m ijT m ki) , 



(7) 



and R is its trace R = 7* i Ru- 
in the BSSN formalism 0,0, the above ADM equa- 
tions are rewritten by introducing the conformally related 
metric 7y 



7y — e 7y') 



(8) 



P = 7i fe f} fe = -7^., 



(12) 



are regarded as independent variables in this formulation. 

The evolution equations of BSSN formulation can be 
written as 



d 
~dt 



1 



= -2aA tj , 
jK = af^A^ + ^K 2 ) - : ' l D l D l n 



(13) 
(14) 



Jt Aj 



.!,, = a (KAij - 2A ik A k j 

+e- 4 *(aR lJ ~D l D 1 a) TF , 



(15) 



d t T l = 2a { f) k A> k - -7 ij X- + QA ij 



1M '•'«.,. + l 3k P\ 3 k + \l ij P k ,ik + P j f \ 3 



(16) 



Here the superscript TF denotes the trace-free part of a 
tensor. The Ricci tensor R^ can be written as a sum of 
two pieces 



Rij — Rij + i2y , 



(17) 



where i?f is given by 



Rfj = -2D l D J ^-2%D k D k <j) 

+Ab l $b 3 4>-A% J b l <t>bi(t )l (is) 

while, with the help of the T l , R^ can be expressed as 

Rij = --ji mn iij,mn+ lk{i^ k j) +f fc f(ij)fc 

+7 m ™ (vf fe m (if j^ kn + t k in f kmj ^j . (19) 
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The new variables are tensor densities, so that their Lie 
derivatives arc 



Cptf, = f3 k cj>, k + ^f3 k , k , 



(20) 



= P%j,k+2j k(i p k !j) --%f3 k , k , (21) 



C K = f3 k K ik , 

ok 



(22) 



C p A l3 = p*Ai j>k + 2A kii /3* tj) - -A %3 (3\ k . (23) 

The Hamiltonian and momentum constraints J5J and (0 
can be rewritten as 

H = (R - 8£) 4 A0 - 8.DVA 

+ \k 2 - A l3 A lj = (24) 
M* = DjA ij + 6A ij <P tj - \f j K d = 0, (25) 
where R = j v Rij. 



III. BLACK HOLES IN KERR-SCHILD 
COORDINATES 

The ingoing Kerr-Schild form of the Kerr metric is 
given bv |2ll39l| 



ds 2 = {j]^ v + 2H£^)dx fl dx v , 



(26) 



where r/^ = diag(— 1, 1, 1, 1) is the Minkowski metric 
in Cartesian coordinates, and H a scalar function. The 
vector is null both with respect to r/^ and , 



ifl^ = g^£^£ v = 0, 



(27) 



and we have £\ = tli. The general Kerr-Schild BH met- 
ric has 



H 



Mr 



and 



rx + ay ry — ax z 
r z + or r z + a z r 



(28) 



(29) 



Here M is the mass of the Kerr BH, a = J/M is the 
specific angular momentum of the BH, and r and 6 are 
auxiliary spheroidal coordinates defined in terms of the 
Cartesian coordinates by 



x 2 + y 2 : - ^ 
r 2 + a 2 r 2 



(30) 



and 2 = rcos9. The event horizon of the BH is located 
at 



r ch = M + 



(31) 



Comparing l|26|) with the ADM metric one identifies 
the lapse function a, shift vector /3, and the spatial 3- 
metric 7y as 



a = (1 + 2//)- 1 / 2 , 



7ij 



2HL 



(32) 
(33) 
(34) 



We can see here that these variables all extend smoothly 
through the horizon and their gradients near the hori- 
zon are well-behaved.. Given these metric quantities, the 
extrinsic curvature Kij can be computed from Q 



Kj 



K 



2aH( k (£ i £ j H <k + 2H£ (il d k £ lj) ) 
+2a(£ (l d n H + Hd {t £ j} ), (35) 
2a 3 (l + H)PH ti + 2aH£\ l . (36) 



In the static case a = 0, the above expressions reduce 
to the Schwarzschild expressions in in-going Eddington- 
Finkelstein form Hfl 



H = 


M/r, 


^ = 




Kij = 


2M 


r 4 (l + 2Af/r)V2 


K = 


2M 


r 2 (l + 2M/r) 3 / 2 



M 



(37) 



(l + 3¥/r). 



x 2 + y 2 + z 2 . 



where M is the total mass-energy and r 

IV. ADJUSTING THE BSSN EQUATIONS 

For a solution of the BSSN equations to be equivalent 
with a solution of the ADM equations, the new auxiliary 
variables have to satisfy new constraint equations. In 
particular, Ay has to be traceless 



A = f^Aij = 0, 



(38) 



the determinant of the conformally related metric 7^ has 
to be unity 



V = dct(7y) -1 = 0, 



(39) 



and the conformal connection functions T l have to satisfy 
the identity l(T2*j) 



f 1 _ fyjkfi 



jk 



0. 



(40) 



These conditions can be viewed as new constraints, in 
addition to the Hamiltonian and momentum constraints 
PH l and P5)l. 

In an unconstrained evolution calculation, the con- 
straints are monitored only as a code check. It may be 
advantageous, however, either to enforce at least some of 
the constraints during the evolution, or to add evolution 
constraint equations to the evolution equations. 
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Alcubierre and Bragmann |29| . for example, found im- 
proved stability properties when the algebraic constraint 
is enforced. This can be achieved by replacing Ay- 
wit h 



AB1: A, 



An — 



ry. .^ mn A 



(41) 



after each time step. They also found improvements, 
at least in octant symmetry, when using the conformal 
connection function T 1 only in places where it is differen- 
tiated, and instead use, by virtue of 1)40(1 , the contraction 
of the Christoffel symbols everywhere else 



AB2: 



(42) 



Other ways of imposing some of the constraints are 
possible (compare |8lj). While the rule AB1 (equation 
I41|l is appealing in that it treats all components of Ay 
identically, we found particularly stable results by dy- 
namically evolving only five of the six components of 
Aij, and computing the zz component using the alge- 
braic constraint 



To suppress exponential growth of T l 1 \ should have the 
same sign as (3 3 .. 

For sufficiently slowly rotating black holes {J/M < 
0.7 M), we were able to follow the evolution for thou- 
sands of M without encountering any instabilities (see 
Section lv"fl and TableHJ). For more rapidly rotating black 
holes (J/M ~ 0.9M) we experimented with adjustments 
suggested by Yoneda and Shinkai |37j , namely 



7y = rhs of (fH|) — KiaHjr 



and 



dt 



^A i3 = rhs of (O - K 2 ae" 40 7i^ fc , fc , 



(46) 



(47) 



where k,\ and Ki are positive numbers. We did find im- 
provements with the adjustment (|46[) . but not with 147|) 
(see Table |D . 



NUMERICAL IMPLEMENTATION 



A z 



A x x + A y y + A xz r z + AyrfV 



(43) 



We similarly determine ^ zz from the other five metric 
components using the algebraic constraint i|39|) 



~ 2 

lyylxz 



^^ixy^fyz^fxz 



' ^fxx^fyz 



'fxx'fyy 



Ixy 



(44) 



In all our simulations, we use these two identities instead 
of rule AB1. 

Combining (021 an( l <E3J wi th rule AB2 leads to ex- 
ponentially growing, unstable modes if no symmetry as- 
sumption is used. While it is not clear what exactly 
causes this instability in the absence of symmetry as- 
sumptions, one possible source is the last term in the 
evolution equation JT^ for f\ (2/3)P/3 j \ y If /3 j j > 0, 
as for example for single black holes in Kerr-Schild coor- 
dinates, then this term may lead to exponential growth 
of any numerical error in T l (compare a similar discus- 
sion in |3l|). The sign of this term can be reversed by 
adding a multiple of the product Q 1 j3 J a to the evolution 
equation i fTfijl . specifically 



d t f l = rhs of DU) 



2a f T) k A* 



2, 



2A i3 a 



x+3 )T^ l M-x? 1 



0\ 



(45) 



Here \ 1S a free parameter, which in principle can even 
be chosen dynamically during an evolution calculation. 



Our numerical implementation follows very closely the 
recipe suggested by Alcubierre and Brugmann j2j|. We 
finite-difference the evolution equations using an iterative 
Crank-Nicholson scheme with two corrector steps 41] . 

We use centered differencing everywhere except for the 
advection terms on the shift (terms involving f3 l di). For 
these terms, a second-order upwind scheme is used along 
the shift direction. 

We adopt "1+log" slicing [UlUlIi 



(48) 



d t a = - aK 



to specify the lapse a. The shift (3 l is determined either 
from the analytic solution, or from the "Gamma-driver" 
condition 



dtP = Xd t T l 



(49) 




), where we choose A = 0.05 in our simulations 



On the outer boundaries of the numerical grid we im- 
pose a radiative boundary condition that is imposed on 
the difference between a given variable and its analytic 
value / — /analytic = u{r — t) / V where u is an outgoing 
wave function. We apply this condition to all fields ex- 
cept P which we leave fixed to their analytic values at 
the boundary. 

We experiment with both cubic and spherical excision 
regions. For cubical excision regions, we adopt the recipe 
suggested by [2!| to copy the time derivative of every field 
at the boundary from its value on a neighboring grid- 
point. For surfaces, we copy from the nearest grid-point 
along the outward normal, and for edges and corners from 
the nearest grid-point along the corresponding diagonal. 

A generalization of this copying algorithm for spher- 
ical excision methods has been suggested in 46]. For 
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TABLE I: Input parameters for selected evolutions. For each evolution we list the lapse condition, the shift condition, the 
symmetry used, the parameter \ in equation 14511 . the parameter Ki in equation (1461 . the excision shape, the time when the 
changes in all the variables reach the level of round-off error (if achieved), and the total run time. 



Case" 


J/M 


Lapse 


Shift 6 


Symmetry 


x d 




excision 


Machine 
accuracy 6 


Run 

time^ 


OO 





Ana 


Ana 


Oct 


AB2 





cube 


No 


< 100M 


01 





1 + log 


Ana 


Oct 


AB2 





cube 


1200M 


> 2000M 


02 





1 + log 


r 


Oct 


AB2 





cube 


1000M 


> 2000M 


03 





1 + log 


Ana 


Oct 








cube 


1800M 


> 3000Af 


04 





1 + log 


r 


Oct 








cube 


1200M 


> 2000M 


05 





1 + log 


Ana 


Oct 








sphere 


2200M 


> 3000M 


06 





1 + log 


Ana 


Oct 


1/3 





cube 


1600M 


> 2000M 


EO 





1 + log 


Ana 


Equ 


AB2 





cube 


No 


1900M 


El 





1 + log 


Ana 


Equ 








cube 


1800M 


> 4000Af 


E2 





1 + log 


r 


Equ 








cube 


1200M 


> 3000M 


E3 





1 + log 


Ana 


Equ 


1/3 





cube 


1600M 


> 3000M 


E4 





1 + log 


Ana 


Equ 


1/3 





sphere 


No 


< 200M 


E5 





1 + log 


Ana 


Equ 


3/4 





sphere 


1400M 


> 2200A/ 


Nl 





1 + log 


Ana 


None 


2/3 





cube 


1500M 


> 3000M 


N2 





1 + log 


r 


None 


2/3 





cube 


1000M 


> 2000M 


N3 





1 + log 


F 


None 


2/3 





cube 


700M 


> lOOOAf 


N4 





1 + log 


r 


None 


2/3 





cube 


500M 


> 1000A/ 


Al 


0.7M 


1 + log 


Ana 


Equ 


1/3 





cube 


3000M 


> 4000M 


A2 


0.7 M 


1 + log 


r 


Equ 


1/3 





cube 


2000M 


> 3500M 


A3 


0.9M 


1 + log 


Ana 


Equ 


1/3 





cube 


No 


> lOOOAf 


A4 


0.9M 


1 + log 


Ana 


Equ 


2/3 





cube 


No 


> 3200Af 


A5 


0.9M 


1 + log 


Ana 


Equ 


2/3 


0.1 


cube 




> 3300A/ 


A6 


0.9M 


1 + log 


Ana 


Equ 


2/3 


0.2 


cube 




> 3400M 


A7 


0.9M 


1 + log 


r 


Equ 


2/3 





cube 




> 4600M 


A8 


0.9M 


1 + log 


r 


Equ 


2/3 


0.1 


cube 


4600M 


> 5700Af 



"Cases N3 and N4 have smaller domains; case N4 uses a smaller 
grid spacing (see text). 

Ana, analytic; T: Gamma driver CEa.l49l. 
c Oct, octant; Equ, equatorial. 

d AB2, adopting rule AB2 CEq.liSl instead of using Eq. 1451 . 
Tor the symbol " — " see Fig. 1101 

^The symbol ">" means the run is terminated at this time but 
could continue. 



a boundary grid-point (i,j,k) of the excised sphere we 
take its nearest neighbors along the coordinate axes away 
from the center of the BH, say (i + 1, j, k), + 1, k) 
and (i, j, k + 1). These three points define a plane, and we 
interpolate to the intersection of this plane with the nor- 
mal on the surface of the excised region. If one of these 
three neighbors is also a boundary grid-point, we project 
the normal into the line defined by the remaining two 
neighbor points, and do the interpolation there; if two 
of the neighbors are inside the excised region we directly 
copy the remaining third point. If all three points are 
inside the excised region we average the three diagonal 
points (i + 1, j + l,fc), (i + l,j, k + l) and (i, j + l,fc+l); 
if that is also not successful we finally copy the point 

(i + l,i + i,fc + l). 

We empirically find improved stability if copying along 
the x or y direction is given higher priority than copy- 
ing along the z direction. This asymmetry is probably 
introduced by (|43|l and (|44|l . We therefore use the grid- 
point (i, j, k + l) ((« + 1, j, k + l) and/or (i, j + 1, k + 1)) 



only when the grid-points (i + l,j,k) and (i,j + l,k) 
((i + 1, j + 1, k)) are not available. 



VI. NUMERICAL RESULTS 

Our simulations are summarized in Table [I] For most 
of these simulations we use computational domains of 
size < x,y,z < 12M for octant symmetry, —12M < 
x,y < 12M and < z < 12M for equatorial symmetry, 
and — 12A/ < x,y,z < 12M for no symmetry, with a grid 
spacing of Ax = 0AM. In order to analyze the effect of 
resolution we also performed the two Cases N3 and N4 
on a smaller domain of — 6M < x,y,z < 6M and used a 
resolution of Ax = 0AM for N3 and Ax = 0.2M for N4. 
An additional simulation with a resolution of Ax = 0.8M 
(not included in Table [IJ was used to establish second 
order convergence of our code in regions not influenced 
by the excision surface. We always use a Courant factor 
of 1/4 so that At — Ax/4. We excise cubes of volume 
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Mass/M 



Constraint residual 



0.96 



0.94 



0.92 



0.9 



0.88 



Inner surface + volume 
Outer-surface 











Ham. 




mom._ 









0.07 
0.06 
0.05 
0.04 
0.03 
0.02 



Angular momentum/M 



rms of Af 




1000 2000 3000 
t/M 



1000 2000 3000 
t/M 



FIG. 1: The monitored quantities as functions of time for Case 
03. The upper-left panel compares different integrals for the 
ADM mass. The solid line is obtained by using Eq. 1A12I ; the 
dashed line is obtained by using rhs of Eq. fKTH . The radius 
of the inner surface is 1.5M and the radius of the outer surface 
is 11. 5M. The lower-left panel compares different integrals 
for the angular momentum. The solid line is obtained by 
using Eq. 1A20L the dashed line is obtained by using rhs of 
Eq. iA"T9l . The upper-right panel shows the L2 norms of the 
Hamiltonian constraint TL (the solid line) and the momentum 
constraint M x (the dashed line). The lower-right panel shows 
a log plot of the root mean square (r.m.s.) of the changes in 
the lapse (the solid line) and the trace of extrinsic curvature 
(the dotted line) between consecutive time steps. 



(1M) 3 in octant symmetry, (2M) 2 x 1M in equatorial 
symmetry, and (2M) 3 without symmetry assumptions, 
or spheres of radius 1M. 

We will call a simulation "stable" if changes in dynam- 
ical variables 4>, ^ij, K, Aij, T l , a and (3 l drop to round- 
off error (of about 10 -16 in double-precision), and remain 
at that level for several hundred M. Reaching round-off 
implies that the numerical solution has settled down to 
the equilibrium solution of the finite-difference equations 
(as opposed to the equilibrium solution of the differen- 
tial equations, which is provided as initial data). In all 
our stable runs, besides monitoring the global quantities 
consisting of the ADM mass, the angular momentum J z , 
the L/2 norms of the Hamiltonian constraint 7i, the mo- 
mentum constraint M x , and the Gamma constraint Q x 
violations, we also monitor the changes in the represen- 
tative variables, i.e., 0, a, K, j X x-> and A xx , until they 
reach and remain at the level of round-off error. Once the 
equilibrium solution with a fixed resolution is achieved, 
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110 



-10 



10" 



10" 



10 



-16 



10" 



* 
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1000 2000 

t/M 
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FIG. 2: The r.m.s. of the change in the trace of extrinsic 
curvature between consecutive time steps as functions of time 
in the static cases with octant symmetry. The lines are labeled 
sequentially from right to left. Case O0 uses a analytic lapse 
condition. Recipe AB2 is used in Cases Ol and 02, while 
the modification 1451 is adopted in Cases 03 — 06. In Case 
05 an excision sphere instead of a cube is used. All of these 
cases are stable and the K's all reach their equilibrium values 
except Case O0. Cases 03 and Case 06 are identical except 
that x — was used in 03 and \ — 1/3 in 06, illustrating 
the effect of the modification 145H in stability. 



we never observe any instability growing up from the 
round-off error at later time. 

In Table Q] we tabulate the time at which any run 
reaches equilibrium, and the time after which the sim- 
ulation is terminated. 

Simulations of non-rotating black holes are labeled "O" 
for octant symmetry, "E" for equatorial symmetry, and 
"N" for no symmetry. 



A. Non-rotating Black Holes 

We first compare several evolutions in octant symme- 
try. Cases OO and Ol are identical except that in OO the 
analytical lapse is used, while in Ol "1+log" slicing H48|) 
is used. OO crashes after a short time, while Ol is stable. 

Cases Ol and 02 confirm that Alcubierre and 
Briigmann's [2^| algorithm (including their rule AB2) 
leads to stable evolution in octant symmetry, both for 
analytic shift and the "Gamma-driver" condition. The 
latter consistently allows the numerical solution to reach 
equilibrium in a shorter time than the former. Simula- 
tions 03 — 06 show that similar results can be obtained 
by replacing AB2 with the inclusion of the constraint in 
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FIG. 3: The the r.m.s. of the change in the trace of extrinsic 
curvature between consecutive time steps as functions of time 
in the static cases with equatorial symmetry. The lines are 
labeled sequentially from right to left. The recipe AB2 is 
used in Cases E0. The modification 1451 is used in Cases 
El — E5 with different values of \- I n Case E4 and E5 an 
excision sphere instead of a cube is used. All of these cases are 
stable except Case E0 and E4. In E0 an instability appears at 
t ~ 700A/ and the code crashes at t ~ 1900M. In E4 the code 
becomes unstable at the beginning and crashes at t ~ 140A4". 



equation l|45l) . In Case 05 an excision sphere instead of 
a cube is used. 

Results for case 03 are presented in Fig.Q] The upper- 
left panel shows two different integrations of the ADM 
masses which are derived in Appendix IA H and illustrated 
in Fig. 1111 The dashed line is computed from a sur- 
face integral at large separation (equation (|A11I) '). while 
the solid line is computed from a volume integral plus a 
surface integral over a small sphere enclosing the black 
hole singularity ( equation IA12|) . We choose a radius of 
i?i = 1.5M for the inner surface and R2 — 11. 5M for 
the outer surface. For R2 — > 00 the two mass integrals 
should agree and should yield the analytic value M of 
the initial data. Our two mass integrals agree to within 
about 5 %. Their difference arises both from the finite 
grid spacing, which generates greater error for surface 
integrals vs. volume integrals, as well as spurious effects 
near the outer boundary, which may more seriously affect 
the outer surface integral. Their deviation from unity is 
a measure of error induced by the proximity of the outer 
boundary to the black hole. 

The lower-left panel in Fig.^shows surface and volume 
integrations of the angular momentum, similar to the 
mass integrations explained above (see Appendix IA 2|l . 
The dashed line is computed from the outer surface in- 
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tegral i|A19(l : the solid line is computed from a combi- 
nation of volume integral and inner surface integral (see 
equation {X20j). For both integrations the angular mo- 
mentum is very close to zero, as it is supposed to be. 

The upper-right panel shows the L2 norms of the 
Hamiltonian constraint TL (solid line) and the momen- 
tum constraint Ai x (dashed line). The lower-right panel 
shows a log plot of the root mean square (r.m.s.) of the 
changes in the lapse a (the solid line) and the trace of 
extrinsic curvature K (the dashed line) between consecu- 
tive time steps. The changes in a and K both decrease as 
exponentially damped oscillations until they reach round- 
off error at about t ~ 1800A/. Using the dynamical shift 
condition (|4Tj|> expedites the damping of these changes 
and helps to stabilize the code (compare, for example, 
cases 01 and 02 or cases 03 and 04). In Fig.|21we com- 
pare the r.m.s. of changes in K for all cases in octant 
symmetry. 

Reconfirming the findings of (2^, we were unable to 
obtain stable evolutions with method AB2 if the symme- 
try is relaxed from octant to equatorial. The result of 
such an evolution (E0) is included in Fig. [31 where we 
plot the r.m.s. changes of K for various different cases in 
equatorial symmetry. For the case E0, the changes again 
drop exponentially until t ~ 700M, but at later times 
they increase exponentially. This exponentially growing 
mode can be extrapolated back to about round-off error 
at / = 0, indicating that the mode is triggered by round- 
off error in the initial data. Using the Gamma-driver shift 
condition instead of the analytic shift leads to similarly 
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N3 are the same except the computational domain of case N3 
is half size of case N2 in length. The change of K of case 
N3 drops faster than in K of case N2 because that light cross 
time in case N3 is shorter than in case N2. The settings of 
Case N3 and N4 are the same except the resolution in case 
N4 is two times higher than in case N3. 



unstable evolutions. However, replacing the method AB2 
with our modification (|45() with \ > 0, we recover stable 
evolutions, in which all changes drop exponentially until 
they reach round-off error. 

In Fig. [21 we compare several different cases in equato- 
rial symmetry, analyzing the effect of the shift condition, 
the value of \ m 1451) . and the shape of the excised re- 
gion. As before, we find that using the Gamma-condition 
instead of the analytic shift leads to a more rapid set- 
tling down to equilibrium, and comparing cases El and 
E3 shows that increasing \ nas a similar effect. Com- 
paring cases E3 and E4, we find that changing the exci- 
sion shape from cubic to spherical with the other settings 
unchanged can destabilize the code, indicating that our 
copying method on an excised sphere leads to larger nu- 
merical error than on an excised cube. Stability can be 
restored by increasing the value of \i as f° r case E5. De- 
tails for the case El are presented in Fig. 

We find that larger values of \ are needed to obtain 
long-term stable evolutions if symmetry assumptions are 
completely removed. In all cases with no symmetry x = 
2/3 is used. We compare the results of these runs in 
Figure [SJ The cases N2 and N3 are identical except that 
that the computational domain of N3 is half as large as 
that of N2. We find that N3 settles down to equilibrium 
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faster than N2, which can be understood in terms of the 
shorter light crossing time. Similarly, cases N3 and N4 
are identical except that N4 has twice the resolution of 
N3, leading to smaller errors and again faster approach 
to equilibrium. 

Fig.[5]shows the result of case N2. In this evolution we 
also locate the apparent horizon every 20 time-steps using 
the 3D finder described in [47|. The apparent horizon 
mass A/ah = (A/lQ-rr) 1 / 2 , computed from its area A, is 
included in the upper- left panel in Fig. and agrees with 
the analytic value to within less than 2 %. 



B. Rotating Black Holes 

We now turn to black holes with non-zero angular 
momentum. Figure [7| shows the results of Case Al 
with a moderate value of the specific angular momen- 
tum a = J/M = 0.7M. Comparing with our results for 
non-rotating black holes we find that the constraint vio- 
lations are larger, and that it takes a longer time for the 
solution to settle down to equilibrium (t ~ 3000M for 
Al). The mass integrals have larger oscillations at early 
times than those for non-rotating black holes, but settle 
down to similar values as before. We also find that the 
values for the angular momentum settle down to values 
within about 2 % of the analytically correct ones. In 
Fig.[S]we compare the r.m.s. of changes in K for Al and 
A2, which again demonstrates that the Gamma-driver 



9 



Mass/M 



Constraint residual 



0.96 
0.95 
0.94 
0.93 
0.92 
0.91 
0.9 
0.89 

0.75 
0.74 
0.73 
0.72 
0.71 
0.7 
0.69 



Inner- 


surface + volume 


Outer 

III. 


-surface 


w~ 










Angular momentum/M 



rms of Af 



Inner- 


surface + volume 


Outer 

1 

-1 


-surface 


1 

h 

''! 

|>k' 




jl 

! 




! 









10- 4 




a 


io- 8 




trK 


io- 8 






io- 10 


\ 




















* 


icr 16 



1000 2000 3000 4000 1000 2000 3000 4000 
t/M t/M 



FIG. 7: The monitored quantities as functions of time for 
Case Al. Labeling is the same as in Fig. 



shift condition helps to stabilize the evolution. 

Instabilities become even harder to control for a = 
0.9M. The larger angular momentum leads to larger 
numerical error, which by itself makes the simulations 
more demanding. In addition, the event horizon for more 
rapidly rotating black holes is smaller. For a = 0.7M, our 
excision cube of physical side length 2M (corresponding 
to a volume of (2M) 2 x 1M in equatorial symmetry), just 
barely fits inside the event horizon. For a = 0.9M, we 
decreased the size of the excision cube to (1.6M) 2 x 0.8M 
so that it does not protrude from the event horizon. 
However, we were unable to obtain stable evolution with 
the reduced size of the excised region, which is probably 
caused by the increasingly large gradients of the gravi- 
tational fields close to the black hole's central singular- 
ity. Interestingly, we were able to obtain stable evolution 
when we left the excision surface at its original size of 
(2A/) 2 x 1M. This surface protrudes from the horizon 
by small amounts at the corners of the cube, which in- 
troduces errors into the solution. However, these errors 
occur very close to the horizon and hardly affect the so- 
lution in the asymptotic region at all. 

In Fig. 1101 we compare the r.m.s. of the change in K 
for different parameter settings. Case A3, in which the 
setting is the same as case Al except increasing a from 
0.7M to 0.9M, shows an exponentially growing mode, 
indicating an instability. The unstable mode still cannot 
be suppressed in case A4 in which^ has been increased 
from 1/3 to 2/3. We find that this instability can be 
controlled by adding the Hamiltonian to the evolution 
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equation for the spatial metric l(lfj|l as suggested by [13 
(cases A5 and A6). Changing from analytical shift to the 
Gamma-driver shift condition turns out to be more effec- 
tive (case A7), and, not surprisingly, we find the fastest 
decay of changes in K by combining both methods (case 
A8). We show details of case A8 in Fig. |5] 



VII. SUMMARY 

We experiment with various modifications of the BSSN 
formulation and study their effect on the stability of nu- 
merical evolution calculations of static and rotating black 
holes. We force the determinant of the conformally re- 
lated metric to be unity and the trace of the tracclcss part 
of the extrinsic curvature to be zero. We modify the evo- 
lution equation for the new auxiliary conformal connec- 
tion functions by adding their constraint equation, and 
also experiment with adjustments of the other evolution 
equations suggested by |37| . 

Most importantly, we find that an instability that 

arises when octant symmetry is relaxed mis eh is si 

can be overcome when the above modifications are em- 
ployed. We demonstrate that both static and moderately 
rapidly rotating black holes can be evolved stably with- 
out encountering any growing modes. Any changes in 
grid functions settle down to round-off error and remain 
there for several 1000 M. 

We find that the dynamically enforced "Gamma- 
driver" spatial gauge condition for the shift leads to 
more stable evolution than using the analytical shift. We 
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also find that cubical excision surfaces, which are more 
straight-forward to implement in Cartesian coordinates, 
work better than spherical excision surfaces. 

While our modifications to not solve all stability prob- 
lems (e.g. for the most extreme rapidly rotating black 
holes), we believe that they lead to significant improve- 
ments that may be a helpful step towards simulations of 
binary black holes and their coalescence. 
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APPENDIX A: EVALUATION OF THE ADM 
MASS AND ANGULAR MOMENTUM 

1. ADM Mass Integration 

The ADM mass is defined in terms of a surface in- 
tegral at spatial infinity. In numerical simulations, this 
integral can be approximated by an integral evaluated on 
a surface near the outer boundaries of the grid (9^2 in 
Fig. Ill|) . To avoid spurious effects of the outer bound- 
aries, it is often desirable to convert this surface integral 
into a volume integral. For black hole spacetimes there is 
the additional complication of singularities in the black 
hole interiors. In this Appendix, we show how a region 
inside the grid can be excluded, so that the mass can 
be computed from a volume integral over the outer re- 
gion and a surface integral over the boundary of excised 
interior regions. 

In Cartesian coordinates, the ADM mass is defined by 
a surface integral at spatial infinity 



run, j 



)dS l 



(Al) 



where dSi = (l/2)j 1 / 2 e i jkdx : >dx k is the surface element 
and £ijk the Levi-Civita alternating symbol. We now 
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FIG. 11: The diagram illustrates the relation between the 
volume integral on the volume Q and the surface integrals on 
the boundaries dQi and dfl2- 



as illustrated in Fig. ^2 For spacetimes with more than 
one black hole, several interior surfaces can be added on 
the right hand side. 

Applying Gauss' law to the first term of the right hand 
side of equation (|A5() yields 



on. 



(f 1 - v .uis, 



[di(r - r»j) + r kl (T k - t Jk 3 )]^d 
(F - r-" ; ;,/.s- 

{R + rt 3 l3 -r 3k t Jlk )^jd 3 x, 



(A7) 



perform a conformal decomposition 

7ij = ^lij ■ 



(A2) 



where we distinguish from the decomposition (jSJ) in order 
to allow for conformally related metrics 7^ with a deter- 
minant different from unity. Assuming the asymptotic 
behavior 



and 



^ ~ 1 + O(-) when r — > 00 (A3) 
r 



lij ~ + when r — > 00, (A4) 



we can rewrite (|A1(I as 

M = r^j> r 2 l lm i 3n [^{l m n, 3 -l 3 ^m) 

)] dSi 

1D7T J qq 

= 17T / ( f 1 " r^)d% - ^ / Styd$. (A5) 

Here the conformal surface element is defined as dSi = 
(1/2)^ 1 I 2 ^ eijkdx 3 dx k , we use the abbreviations f l = 
^ k r l jk and f tJ k = 7 j7 f; fc , and L>i is the three-covariant 
derivative with respect to the metric 7y . 

Using Gauss' law, the surface integral l|A5|) can be con- 
verted into a volume integral. For spacetimes containing 
a black hole, we can exclude an interior region from the 
volume integration and write the integral over an outer 
surface dQ,2 as a sum of a volume integral over f2 and a 
surface integral over an inner surface dfli, 

[ u l dSi= [ d i (^u i )d 3 x + [ u l dSi (A6) 



Similarly, the second term on the right hand side of equa- 
tion i|A5J) yields 



on. 



D 2 ip^d 3 x+f DtydSi, (A8) 



where D 2 = j l3 DiDj. Collecting terms, the ADM 
mass (|A5(I can now be written 

M = — |- <f (r - fi'j - SDtydSi 
16tt J dn2 



16tt 



(A9) 



— <f> (T* - T^j - 8D l ^)dSi 



We now adopt the decomposition of the BSSN formal- 
ism (Section[nJ), in which the conformal factor is written 
as ip = e ^ an d the conformally related metric 7^ = 7y 
is assumed to have determinant 7 = 1 so that T 3 k - = 0. 
Together with the Hamiltonian constraint (|24|) 



£> V = —R + —K 2 - - — AijA 13 - 27re 5 V (A10) 
8 12 8 1 

where we have included the mass-energy source p = 
n il n v T lLV for completeness, the ADM mass <|A9|1 becomes 



M = — |— <f> (F - S&e^dSi 



(All) 



16tt 



d 6 x 



e 5 *( 167rp + AijA 13 - -K' 



16tt 



-r iife r iife 

(F - S&e^dSi, 



(1 - e*)JZ 



(A12) 



12 



where dSi = (l/2)eijkdx-'dx since 7=1. 

We also note that if the conformally related metric falls 
off faster than inversely with r 

% = % + 0(l/r l+a ), a > (A13) 

(compare |48|). the first term in equation l|A5(l vanishes 
and the mass integral reduces to 



M 



1 

27r Jdn 
1 

16tt 



d 3 a 



10-/; • ,1, - -K' 



1 



-—& D l e*dSi, 

2 77 



9fii 



-e*R 
(A14) 



Since the Kerr-Schild metric (|2(j[l does not satisfy the 
fall-off condition l|A13|) , equation (|A12(1 has to be used to 
evaluate its mass (equation (I A 14(1 would yield the incor- 
rect result M/3 for the static Kerr-Schild metric). We 
have also found empirically that, even for metrics for 
which (|A14p is appropriate, (| Al 2|) yields more accurate 
values for the mass than (|A14J| in dynamical evolution 
calculations 1491. 



= e y - fc / xPe^A l k dSi + ey* f {x> e 6 * A 1 k ) jd 3 x . 

(A16) 



The volume integral in Eq. I|A16(I is 



{x'e^A l k )^d 3 x 



5he^A l k +x>{e^A l k )j 



d 3 i 



(A17) 



^A j k + xWtie^A 1 ^ + e^x j t n H A l 



d 3 x 
d 3 x 



e^Ai k + xW^A'k) - -e 6 ^4f", t 



where we have used 7=1. With the momentum con- 
straint, 



D^&i) = (~DiK + 8tvsA , (A18) 

where we have included the momentum density Si = 
7i A1 n w T /il/ for completeness, the volume integral (|A18|) 
can be rewritten 



2. ADM Angular Momentum Integration 

We define the angular momentum J 1 as 

Ji=^-e ij k I ^A l k dS t . (A15) 

(compare 0|), where the indices of ey & are raised 
and lowered with the flat metric 8ij . Since the integrand 
is evaluated at r — > oo, we can replace dSi = e^dSi and 
A l j — A l j and use Gauss' law to obtain 



87TJ, 



x?e b +A l k dSi 



Ji 



1 

8^' 



8tt v 



8tt 



(A19) 



an 2 



e 60 l^Vk + -x j D k K 
~x> : A tn d k j ln + 87rx J s k 



rf ;> 



e^x^A'kdS,. 



(A20) 



on. 



This expression is equivalent to Eq. (2.25) of ||. Note 
that the first term does not vanish identically, since in- 
dices of Cijn are raised with §ij while indices of A> n are 
lowered with 7y . 
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